Introduction to Raster Data

Raster Basics

# Load digital surface model for Harvard Forest
# 1m GEOTIF
# UTM Zone 18
# can navigatre through subfolders by hitting tab
HARV_dsmCrop_info <- GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")

HARV_dsmCrop_info
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##     Bmin   Bmax    Bmean      Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area
# Load the GEOTIF as a raster
DSM_HARV <- raster(here::here("data", "NEON-DS-Airborne-Remote-Sensing", "NEON-DS-Airborne-Remote-Sensing", "HARV", "DSM", "HARV_dsmCrop.tif"))

DSM_HARV
## class      : RasterLayer 
## dimensions : 1367, 1697, 2319799  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/User/Documents/R/spatial-in-r_2019-11/data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif 
## names      : HARV_dsmCrop 
## values     : 305.07, 416.07  (min, max)
# call summary of raster
summary(DSM_HARV) # using 4.31% of cells
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (4.31% of all cells)
##         HARV_dsmCrop
## Min.          305.32
## 1st Qu.       345.49
## Median        359.48
## 3rd Qu.       374.12
## Max.          414.43
## NA's            0.00
summary(DSM_HARV, maxsamp=ncell(DSM_HARV)) # using all the cells
##         HARV_dsmCrop
## Min.          305.07
## 1st Qu.       345.59
## Median        359.67
## 3rd Qu.       374.28
## Max.          416.07
## NA's            0.00
#DSM_HARV is currently a raster... we want a data frame to plot in ggplot

DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE) %>%
  mutate(z = HARV_dsmCrop)
# this df has more than 2 million rows
# consider whether a raster really needs to be a data frame...

ggplot() +
  geom_raster(data = DSM_HARV_df, aes(x = x, y = y, fill = z)) +
  scale_fill_viridis_c() +
  coord_quickmap()

plot(DSM_HARV) # makes a gross looking plot of the raster, but it's paster than ggplot for this purpose

Coordinate Reference Systems

Challenge

Use the output from the GDALinfo() function to find out what NoDataValue is used for our DSM_HARV dataset.

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##     Bmin   Bmax    Bmean      Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area

Found in first data.frame output from GDALinfo() NoDataValue = -9999

Histogram of raster

ggplot() +
    geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Plot Raster Data in R

Reproject Raster Data in R

Raster Projection in R

DTM_HARV <- raster(here::here("data", "NEON-DS-Airborne-Remote-Sensing", "NEON-DS-Airborne-Remote-Sensing", "HARV", "DTM", "HARV_dtmCrop.tif")) # WGS84 using UTM Zone 18

DTM_HARV
## class      : RasterLayer 
## dimensions : 1367, 1697, 2319799  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/User/Documents/R/spatial-in-r_2019-11/data/NEON-DS-Airborne-Remote-Sensing/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif 
## names      : HARV_dtmCrop 
## values     : 304.56, 389.82  (min, max)
DTM_hill_HARV <- raster(here::here("data", "NEON-DS-Airborne-Remote-Sensing", "NEON-DS-Airborne-Remote-Sensing", "HARV", "DTM", "HARV_DTMhill_WGS84.tif")) # WGS84 using lat-long

crs(DTM_hill_HARV)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Reprojecting the hillshade file

Reprojecting changes the data!

DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV, crs = crs(DTM_HARV))
crs(DTM_hill_UTMZ18N_HARV)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
extent(DTM_hill_UTMZ18N_HARV)
## class      : Extent 
## xmin       : 731397.3 
## xmax       : 733205.3 
## ymin       : 4712403 
## ymax       : 4713907
crs(DTM_HARV)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

Math with Rasters: Canopy Height Model

Creating a canopy height model by subtracting digital terrain model from digital surface model

# double check that the projections are the same
crs(DTM_HARV)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
crs(DSM_HARV)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# check that the extents are the same
extent(DTM_HARV)
## class      : Extent 
## xmin       : 731453 
## xmax       : 733150 
## ymin       : 4712471 
## ymax       : 4713838
extent(DSM_HARV)
## class      : Extent 
## xmin       : 731453 
## xmax       : 733150 
## ymin       : 4712471 
## ymax       : 4713838
# they are
# if they weren't, we'd use clip()

# and check them out visually
plot(DTM_HARV) # smoother and lower elevation

plot(DSM_HARV) # rougher and higher elevation

#### With direct subtraction

# quick and dirty canopy height model
CHM <-  DSM_HARV - DTM_HARV
plot(CHM)

#### With overlay() overlay() takes two rasters and a function telling it what to do with them

CHM_ov_HARV <- overlay(DSM_HARV, DTM_HARV, fun = function(r1, r2) {return(r1 - r2)})

plot(CHM_ov_HARV)